C~ V 


# 0*7 


NASA Technical Memorandum 110340 



High Accuracy Evaluation of the Finite 
Fourier Transform Using Sampled Data 


Eugene A. Morelli 

Langley Research Center, Hampton, Virginia 


June 1997 


National Aeronautics and 
Space Administration 
Langley Research Center 
Hampton, Virginia 23681-0001 




Table of Contents 


Abstract iii 

Nomenclature iv 

I. Introduction 1 

II. Theoretical Development 2 

III. Examples 9 

IV. Summary 1 1 

V. Acknowledgments 11 

VI. References 12 

Appendix - Derivation of the Cubic Interpolation Formulae 13 

Figures 19 




Abstract 


Many system identification and signal processing procedures can be done advantageously in 
the frequency domain. A required preliminary step for this approach is the transformation of 
sampled time domain data into the frequency domain. The analytical tool used for this 
transformation is the finite Fourier transform. Inaccuracy in the transformation can degrade system 
identification and signal processing results. This work presents a method for evaluating the finite 
Fourier transform using cubic interpolation of sampled time domain data for high accuracy, and the 
chirp z-transform for arbitrary frequency resolution. The accuracy of the technique is 
demonstrated in example cases where the transformation can be evaluated analytically. Arbitrary 
frequency resolution is shown to be important for capturing details of the data in the frequency 
domain. The technique is demonstrated using flight test data from a longitudinal maneuver of the 
F-18 High Alpha Research Vehicle. 
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Nomenclature 


/ cyclic frequency, Hz 

FFT Fast Fourier Transform 

i time index 

j imaginary number = V-l 

k frequency index 

M number of discrete frequencies 

N number of sampling intervals 

t time, sec 

T data record length, sec 

P n nth order Lagrange interpolating polynomial 

x{t) real function of time t 

Xj sampled x(t) at time iAt 

x k discrete Fourier transform of jc, for frequency f k 

x(f) finite Fourier transform of x(t) for frequency/ 
W interior point cubic weighting function 

a angle of attack, deg 

At sampling interval, sec 

(p endpoint interpolation function 

y endpoint cubic weighting function 

6 discrete Fourier transform angular step 

co frequency, rad / sec 

y/ interior point interpolation function 

superscripts 

Fourier transform 
* complex conjugate 


iv 



I. Introduction 


Data analysis and linear parametric modeling can be carried out in either the time domain or the 
frequency domain 1 . Frequency domain analysis has many advantages, including physical insight, 
direct applicability to control system design, and lower dimensionality for model parameter 
estimation, among others 1-3 . Frequency domain techniques based on spectral estimation 3-5 , output 
error 1-3 , and equation error 1 ’ 6 ’ 7 have been used successfully for flight test data analysis and 
modeling. The basis for all of these frequency domain methods is the finite Fourier transform, 
which is the mechanism for transforming time domain data to the frequency domain. Any errors in 
the transformation from time to frequency domain affect the accuracy of the raw data in the 
frequency domain, which in turn impacts data analysis and modeling results. 

Typically, transformation from the time domain to frequency domain is brought about using 
sampled time domain data and a simple Euler approximation of the finite Fourier integral. This 
calculation, called the discrete Fourier transform, can be implemented with a Fast Fourier 
Transform (FFT) algorithm 8 for increased computational speed. The simple Euler approximation 
makes the discrete Fourier transform increasingly inaccurate with increasing frequency and 
sampling time. Accuracy of transformation to the frequency domain can be improved using 
various quadrature methods 8 , which involve evaluating the Fourier integral separately for each 
individual frequency of interest. 

For a given data record length, the discrete Fourier transform constrains the frequency 
resolution to equal the reciprocal of the record length. As the record length decreases, frequency 
resolution can become coarse, leading to omission of important data features in the frequency 
domain. One solution to this problem is to pad the time domain data with zeros and thus artificially 
increase the data record length, but this is inconvenient and increases processing time. It is also 
possible to interpolate values in the frequency domain, but this is also inconvenient and can be 
inaccurate as well. When the frequency band of interest is relatively small, most of the processing 
involved in the discrete Fourier transform is for frequencies outside the range of interest, since the 
computations include all frequencies up to the Nyquist frequency. 

The purpose of this work is to present a convenient and efficient method for computing the 
finite Fourier transform with high accuracy using sampled data, while allowing a wide range of 
choice for the frequency band and frequency resolution. This capability gives a solid underpinning 
to the many frequency domain techniques developed and renders those methods commensurately 
more accurate and flexible. 

The next section outlines the theory for the flexible, high accuracy finite Fourier transform 
method. Following that, some realistic simulation examples are presented to demonstrate some of 
the issues discussed in the theory. Finally, an example using measured flight test data from the 
F-18 High Alpha Research Vehicle is included. 
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II. Theoretical Development 

The Fourier transform of a continuous scalar time function x(t) on a finite time interval [0, T ] 
is called the finite Fourier transform, and is defined by 


*(/) = £ x(t)e- J2 * f 'dt 
or 

x(co) = j Q r x( t ) e~i m dt 

where 

co = 2 k f 


( 1 ) 

( 2 ) 

( 3 ) 


When the time function x(t) is sampled at discrete evenly spaced time intervals, the finite 
Fourier transform can be approximated by 


N - 1 

*(/) = A f £ x, e ~ j2lcft ‘ (4) 

[=0 

where 

t t =iAt , x, , = x(tj) = x(iAt) i =0,1,2,..., N (5) 

and 

A t = — (6) 

N K ’ 


Eq. (4) is a simple Euler approximation for the integral in Eq. (1), using the first N discrete 
values of the continuous time function x(t). Note from Eq. (5) that the total number of samples is 
N + l. The discrete Fourier transform can be used to evaluate the summation in Eq. (4) at specific 
frequencies. The discrete Fourier transform is defined by 


N - 1 

x k a V x,e~ J2 * f & k =0,1,2, (7) 

i=0 


2 



where 


and 


k _ k 
NAt ~ T 


k =0,1,2, M-\ (8) 


r n /2 

M = \ 

l(* + l )/2 


for N even 
for N odd 


(9) 


The values given for M in Eq. (9) represent the fundamental limitation that frequencies 
contained in a sampled time history must fall within the frequency band [0,/ n ) (i.e., below /„), 
where /„ is the Nyquist frequency, defined as one half the sampling frequency: 


/«- 


1 

2 At 


For the discrete frequencies f k in Eq. (8), using Eq. (7) in Eq. (4) gives 


( 10 ) 


x(f k )~Atx k k =0,1,2, ...,M-1 (11) 


The discrete Fourier transform in Eq. (7) can be expressed in another form by substituting 
from Eqs. (5) and (8) to obtain 


x k = X x i e jl7Cik/N 

i = o 


it = 0 , 1 , 2 ,..., A /-1 (12) 


Now define 

9 k 

6 = co k At = 2nf k At = — k =0,1,2,..., M-l (13) 

N 


where the dependence of 0on k is not shown explicitly by the notation. Using Eq. (13), Eq. (12) 
can be written as 
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k =0,1,2,..., M-l (14) 


JV-l 

*k = X x « e ' y0 ' 


1=0 


Combining Eqs. (12)-(14), and (4), 


x(0) = x(2nf k At) = Atx k k =0,1,2, ... , M -l (15) 


When N is an integer power of 2, FFT algorithms 8 can be used to efficiently carry out the 
summation in Eq. (12). 

Eqs. (7) and (8) indicate that the discrete Fourier transform is defined for frequencies which 
correspond to an integral number of cycles for each time period T. Since the finite Fourier 
transform is an integration over the time interval [0, T], the values computed using the 
approximation in Eq. (15) can be small, and thus susceptible to round-off error accumulated 
during the summation in Eq. (12). In addition, since Eq. (15) represents a simple Euler 
approximation to the integral in Eq. ( 1 ), high frequencies f k or large sample times At will 

introduce additional inaccuracy because of the oscillatory nature of the integrand. The latter 
difficulty was recognized as early as 1928 by Filon 9 , who proposed an interpolation scheme for 
the sampled time domain data to improve the accuracy of finite Fourier integrals. The same basic 
approach is taken here, using a cubic Lagrange polynomial interpolation due to Press et al. 8 . The 
resulting calculation uses the discrete Fourier transform with weighting factors as follows: 


x(0) = At{ W(6)[x k +x N e-J eN ] 

+ 7o(0)*o + 7i(0)*i + 72(0)*2 + Ys(0)x 3 (16) 

+ e j9TIAl '[7o(0)*jv + 7i*(0)*am + Yi{ 0 )x N _ 2 + 73*(0)*n-3]} 


where 


W(0) = 


6 + fl 2 
30 4 J 


-4cos6 + cos20) = 1 


720 15120 


0 6 


(17) 
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and 


(-42 + 50 2 ) + (6 + e 2 )(8cos9-coi2e) (-l20 + 60 3 ) + (6 + 0’-)j<n20 

ToW - j — a 


.-? + +* 2 + J5L e 4_. '69 


3 45 15120 226800 


e 6 -jd(— +— e 2 — — 0 4 + 
V45 105 2835 


( 18 ) 


86 V 


467775 J 


14(3- 0 2 ) — 7(6 + 6 2 }cos 0 300-5(6 + 0 2S )sinO 

Yi(0)- 604 j 60^ 


7 7 0 2 +-J-0 4 7 


24 180 3456 259200 


0 6 -jd ( — — 0 

l 72 168 


(19) 


2 + — 0 4 - ^ 06 
72576 5987520 


-4(3- 0 2 ) + 2(6 + 0 2 }cos 0 - 120 + 2(6 + 0 2 )«n0 

YiX )- 3 0 4 ~ J 304 


i + +0 2 . 

6 45 6048 


— 0 4 + — - — 0 6 — jd (— — H — — 0 2 - 
)48 64800 l 90 210 


( 20 ) 


11 ■ 0 4 + — 0 ^ 


90720 7484400 


2(3 - 0 2 ) - (6 + 0 2 )co 5 0 60 -(6 + 0 2 )«'/i0 

604 7 604 


(21) 


= — — 0 2 + — — — 0 4 — 


24 180 24192 259200 


- — 0 6 - jo( 
)200 V 


7 1 02+-^L_0 4 - 13 06 


360 840 362880 29937600 


Detailed derivation of the expressions in Eqs. (16)-(21) can be found in the Appendix. The 
series expansions for W(0), y o (0), 7 i( 0), 7 2 (0), and 7 3 (0) given above are used for small 0, 
where 0 is considered small when its magnitude is less than the largest value of 0 for which 
identical results are obtained to machine precision with both the analytic expression and the series 
expansion. The series expansions are necessary because of high order cancellations that make the 
analytic expressions inaccurate when 0 is small. 

For many practical problems, such as flight test data analysis, only a relatively small frequency 
band is of interest. It is advantageous to select the frequencies in the finite Fourier transform to be 
closely spaced within the frequency band of interest to ensure that details of the frequency 
spectrum are accurately captured. Arbitrary frequency resolution in a selected frequency band can 
be achieved with the chirp z-transform 10 . The flexibility of the chirp z-transform allows high 
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resolution calculation of the Fourier integral in important regions of the frequency domain, with up 
to N discrete frequencies evenly spaced in a selected frequency band. 

Suppose that M discrete frequencies are selected in the frequency band [/ 0 , f \ ) , 

f k =f 0 +kAf k =0,1,2,..., M -1 (22) 

where 


Af= ih ~/o) 
M 


(23) 


Then Eq. (7) can be written as 


N - 1 N - 1 

x k = Y i x i e~ j2nfk, < = £ x, e~ j27lf o iAt e -j2nkAfiAt k =0,1,2, , M - 1 (24) 
1=0 i =0 


Defining 


<j) 0 = 2Kf 0 At , A<f> = 2nAf At 


(25) 


A = ( z = <?' A * 


(26) 


and combining Eqs. (24)-(26), 


N - 1 N-l 

** = £*,a-'z-‘' = £x, [az‘J 

j =0 i =0 


(27) 


As A: increases, the quantity AZ* represents a contour along the unit circle in the z plane. The 
quantity <p 0 is the angular location on the unit circle associated with the starting frequency / 0 , and 
A<j> represents the incremental angular step along the unit circle in the z plane for each frequency 
increment A/, corresponding to each increment in k. Figure 1 illustrates the idea. Comparing 
Eqs. (14) and (24), and using Eq. (25), 


9=<l> 0 + kA<t> = 2nAt(f 0 + kAf) k =0,1,2,..., M -1 (28) 
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Suppose the initial frequency is chosen as zero, 


fo —Q > 0o - 0 > A — l 


(29) 


and the frequency increment is chosen as 



1 

NAt 



, Z = e j lK ! N 


(30) 


Then if the number of discrete frequencies M is chosen from Eq. (9), the chirp z-transform is 
identical to the discrete Fourier transform of Eq. (12). Using Eqs. (29) and (30) in Eq. (28), 


6 = kA<p 


Ink 

N 


k =0,1,2,..., Af-1 (31) 


The choices in Eqs. (29) and (30) correspond to the case where the contour in the z plane 
consists of equal angular increments around the entire unit circle, which is identical to that for the 
discrete Fourier transform of Eq. (12), cf. Eqs. (13) and (31). Eqs. (7)-(9) indicate that the 
contour in the z plane for the discrete Fourier transform completes nearly a half circle from 0 = 0 
up to but not including 0 = n, which corresponds to the Nyquist frequency. For the chirp 
z-transform, the starting frequency and frequency increment are not restricted to the values given in 
Eqs. (29)-(30), but can be chosen arbitrarily. This translates to small angular increments over an 
arc of the unit circle in the z plane. The chirp z-transform is also called the zoom Fourier 
transform, because of the capability to achieve very fine resolution in the frequency domain for a 
selected frequency range. 

A disadvantage of using the frequencies associated with the discrete Fourier transform 
(specified in Eqs. (8) and (9)) is that the spacing of these discrete frequencies depends on the data 
record length T. For short data records, the frequency resolution can be coarse, leading to a loss 
of fine detail in the frequency domain, with consequent degraded data analysis and modeling 
results. In extreme cases, it is possible for the discrete frequencies to be spaced such that 
significant features in the frequency domain are missed. In addition, the frequency points used in 
the discrete Fourier transform are evenly spaced between zero frequency and the Nyquist 
frequency f n . For a small bandwidth of interest, this means that much of the computation is for 

frequencies outside the range of interest. The resulting number of usable data points in the 
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frequency domain can therefore be small, again to the detriment of the data analysis and modeling 
results. 

The chirp ^-transform decouples the frequency resolution from the length of the time record 
and can place all calculated frequency points within the frequency band of interest. The cubic 
interpolation corrections given above in Eqs. (16)-(21) apply without modification to the chirp 
z-transform because the chirp z-transform is just a discrete Fourier transform using discrete 
frequencies that can be adjusted arbitrarily, independent of the time record length. The penalty for 
this flexibility is more computation time, but the magnitude of this extra computation time 
continues to shrink as computing machinery becomes faster. It is possible to compute the chirp 
z-transform efficiently for arbitrary N using FFT algorithms 10 . 

In summary, the procedure for computing the high accuracy chirp z cubic approximation to the 
finite Fourier transform is: 


1 . Choose the initial frequency / 0 and the number of frequency domain points M. 

2. Choose the frequency resolution A/, or choose the upper limit of the frequency band f x 
and compute A / from Eq. (23). 

3 . Compute the values of 6 from Eq. (28). 

4. Use the chirp z-transform to carry out the summation in Eq. (14), which is equivalent to 
Eq. (27). 

5 . Compute the weighting factor and endpoint corrections from Eqs. (17)-(21) for each Q. 

6. Use Eq. (16) to compute the high accuracy chirp z cubic approximation to the finite 
Fourier transform for each 6. 
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III. Examples 

Figure 2 is a plot of an example time function, x(t) = 1 - e~ 2t , At = 0. 05, T = 5, which might 
represent the step response of a dynamical system. Figure 3 shows the magnitude and phase of 
the finite Fourier transform of x(r) calculated using the Euler approximation of Eq. (15) and the 

cubic interpolation approximation of Eq. (16), compared to values calculated analytically using 
Eq. (1) for the same frequencies f k from Eq. (8). The analytic values and the cubic 

approximation are very nearly identical, so that the associated traces are indistinguishable in 
Figure 3. The Euler approximation shows some inaccuracy in magnitude, with phase angle 
inaccuracy that increases with frequency. Figure 4 shows the differences in magnitude and phase 
between the analytic values and the approximate values, and Figure 5 is a plot of the percent errors 
based on the analytic values. All results in this work were computed using double precision 
arithmetic. Accuracy of the Euler approximation would be worse using single precision arithmetic 
because of the susceptibility of the calculation to round-off error, discussed previously. Figures 4 
and 5 quantify the errors associated with using the Euler approximation and clearly show that the 
cubic interpolation scheme of Eq. (16) completely corrects the errors. This example demonstrates 
the high accuracy associated with the cubic interpolation scheme compared to the conventional 
Euler approximation. 

Figure 6 is a different example time function, x(t) = 5e~‘ sin(nt), At — 0. 05, T = 5, which 
might represent the impulse response of a dynamical system. The frequency spectrum for this time 
function is mainly associated with the frequency of the sine function, in contrast to the last example 
where the frequency content was dominated by low frequencies near zero. In Figure 7, the Euler 
approximation of Eq. (15) is compared with the cubic interpolation approximation of Eq. (16) 
using the chirp ^-transform to achieve frequency resolution ten times finer than that computed from 
Eq. (8) and used in the Euler approximation. The magnitude plot in Figure 7 shows that the Euler 
approximation does not capture the full detail of the frequency spectrum due to the relatively coarse 
resolution in the frequency domain resulting from Eq. (8). The horizontal scale covers only 0 to 
2 Hz and does not extend to the full Nyquist frequency (10 Hz), in order to illustrate the frequency 
resolution. Figure 8 shows the differences in magnitude and phase between the analytic values and 
the approximations. The comparison of the Euler approximation with analytic values used the 
frequencies from Eq. (8), while the comparison of the chirp z cubic approximation to analytic 
values used the tenfold higher frequency resolution. Figure 8 is therefore an accuracy check 
analogous to Figure 4, and shows that the chirp z cubic approximation is also accurate for the 
higher frequency resolution. At higher frequencies not shown in Figure 8, the phase angle error 
for the Euler approximation increases in a fashion similar to that seen in Figure 4 for the last 
example. This example demonstrates that the chirp z cubic approximation can not only provide 
details of the frequency spectrum that are missed with the coarse frequency resolution of the Euler 
approximation, but also produces values that are highly accurate for the finer frequency resolution 
of the chirp z cubic approximation. 
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Figure 9 shows angle of attack measurements taken during a 15 second longitudinal flight test 
maneuver flown on the F-18 High Alpha Research Vehicle. The s amplin g rate was 50 Hz. 

Figure 10 shows the corresponding magnitude and phase angle results computed using the Euler 
approximation of Eq. (15) and the chirp z cubic approximation of Eq. (16) with a frequency 
resolution ten times finer than that computed from Eq. (8) and used in the Euler approximation. 
The frequency scale covers 0 to 1.4 Hz, because all the rigid body dynamics he in this frequency 
band. In general, it is not necessary to choose the frequency resolution to be an integer multiple of 
the f k from Eq. (8) for the chirp z cubic approximation. In fact, Eqs. (22) and (23) indicate that 

the frequency resolution is arbitrary, with the only limitation being that the frequency steps 
associated with Z in Eq. (26) will be the same for each increment in k. It is also possible to move 
the frequency band of interest anywhere below the Nyquist frequency by choosing / 0 and /j . 

This might be desired to investigate high frequency structural modes, for example. For the chirp z 
cubic approximation, it was therefore possible to place all the computed frequency domain points 
in the frequency band given above, with no wasted calculations. In contrast, the Euler 
approximation computes values for frequencies f k out to the Nyquist frequency (25 Hz), but the 
only values of interest are associated with f k in the frequency band 0 to 1.4 Hz. As in the 
previous example, the chirp z cubic approximation captures details of the frequency spectrum that 
were missed using the coarse frequency resolution of Eq. (8). In signal processing and system 
identification applications, this corresponds to an increase in the number of data points in the 
frequency domain, which generally improves results. In Figure 11, the differences between 
magnitude and phase computed using the Euler approximation and the chirp z cubic approximation 
are plotted for the frequencies used in the Euler approximation ( f k from Eq. (8)). Figure 12 
shows the magnitude and phase differences between the two approximations for all f k from 
Eq. (8) out to the Nyquist frequency. Although it is not possible to assess the accuracy of the 
approximations for flight test data because the analytic values are unknown, it can be seen that the 
phase differences in Figure 12 show an increase with frequency that is s imilar to what was seen in 
the simulation example cases. 
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IV. Summary 

A high accuracy method for computing the finite Fourier transform for sampled time domain 
data was presented and demonstrated with examples. The technique uses cubic interpolation of the 
sampled time data in the integrand of the finite Fourier transform, and employs the chirp 
z-transform to achieve arbitrary resolution in the frequency domain regardless of the length of the 
time record. Accuracy improvement was demonstrated using examples, and the effect of full 
control over frequency resolution independent of the length of the time record was demonstrated. 
Although other quadrature methods can be used to accurately compute the finite Fourier transform, 
the technique presented here has the advantages of convenient and efficient calculation, while at the 
same time achieving high accuracy and arbitrary frequency resolution. 
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A ppendix - Derivation of the Cubic Interpolation Formulae 

To accurately evaluate the finite Fourier transform in Eq. (1), Lagrange polynomials are used to 
interpolate samples of the function x(t ) in the integrand. This avoids the inaccuracy associated 
with the Euler approximation of the integral represented by the calculations in Eqs. (14) and (15). 
Lagrange interpolating polynomials have the property that they pass through each data point. The 
nth order Lagrange interpolating polynomial is given by 8 


where 


p„(<)=£ 

m = 0 


LM r 
<-(«.) ’ 




j = 0 

j*m 


(A-l) 


(A-2) 


For nth order interpolation, n + 1 function values are required and the nth order Lagrange 
interpolating polynomial includes n + 1 nth order functions of t weighted by neighboring discrete 
function values. Cubic interpolation therefore requires four function values, two on either side of 
the interval for which the interpolation is being done. The functions £ m (t) are simple nth order 
polynomials in t. In Figure A-l, the sampled function value jc, at a typical interior point t, is 
shown to be involved in the interpolation when t falls within any of four different intervals, 
namely, (*,•-!»*,•)> (t,, t, +[ ), and (f l+1 ,f l+2 ). The endpoints would be involved in fewer 

intervals. The first and last interval have to use three points on one side (toward the interior of the 
data record) and the endpoint. Outside of these intervals for t, the ith function value does not 
contribute to the interpolation. A general interpolation expression for the function x(t) can then be 
written as 


N 

¥(*)+ 

1=0 i=endpoints 


(A-3) 


where the <Pi(t) multiplying the endpoint function values must account for the endpoints being 
included in the first summation as if they were regular interior points. Applying the finite Fourier 
transform to both sides of Eq. (A-3), 
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x«o)~j o 



N 1 

Jo ‘ 


Z, x iVi (O’ 


l' = o 

/=endpoints 


>e 


i m dt 


(A-4) 


Introducing the change of variable 


s = (t-t i )/At 


(A-5) 


for integral of the first summation, so that s t _ j = (/,•_ l - r f )/Af = — 1, s ( = (r, - t { )j At = 0, 

■Sj+i = (t i+ \ - h)/ At = 1 , etc. This is equivalent to mapping t K to zero on the s axis and normalizing 

the sampling interval to unity. For the integral of the second summation, introduce the change of 
variable 


s = t/At 


(A-6) 


which has similar effects. Making the substitutions and using Eq. (5), Eq. (A-4) then becomes 


x(co)~Atj • 

' N 

Y x; 



i=0 

/=endpoints 


(A-7) 


The limits of the integration in s can be set to ±°° because both yr(s) and <p(s) are zero outside a 
finite range of s values. Switching the order of integration and differentiation, and denoting 
6 = co At as in Eq. (13), Eq. (A-7) becomes 


or 


x(co) ~ At 


N 


Wie^e-J" 

1=0 


+ r,(*) - 

i=endpoints 


(A-8) 


x(co) « At < 

3 

x N e j° N + e j Ql 




i=0 

/= endpoints 
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where 


W(0)= P \f/(s)e-j 6s ds 

J—oo 

Yi(0)=\ <Pi( s ) e ~ j0s ds 

J — oo 


(A- 10) 


(A-ll) 


The expressions in Eqs. (A- 10) and (A-ll) need only be evaluated once for a particular 
interpolation scheme. The first summation in Eq. (A-9) is recognized as the discrete Fourier 
transform from Eq. (14). For cubic interpolation, using Eq. (14) in Eq. (A-9) gives 


x(6)~m[ W{e)[x k +x N e~i 6N ] 

+ 7o(0)*o + Yi(0) x i + Y 2( e ) x 2 + 73( 0 )^3 (A- 12) 

+ e j6T /A ‘ [ro(0)x N + yl{0)x N _ 1 + 72 (^)^ v -2 + 73(^)^- 3 ]} 

which is identical to Eq. (16). 

From Eqs. (A-l), (A-2), and (A-5), the Lagrange polynomial in s multiplying an interior point 
x j for the interpolation labeled 1 in Figure A-l is 


y lW= (« + 3X»+ 2X^ + 1) 

6 


Similarly, for the interpolations labeled 2, 3, and 4 in Figure A-l, the Lagrange polynomials in s 
multiplying an interior point x t are 


v , 2 ( j )=(* +2)(i+1)(s - 1) 

(A- 14) 

^( J )= (s+1)(s “ 1)(s “ 2) 

(A- 15) 

—6 

(A- 16) 
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Since an interior point x t is involved in only the four interpolations shown in Figure A-l, 


\f/(s) = ( 5 ) + V 2 (s) + V&) + V 4 (s) 


(A- 17) 


Using Eq. (A- 17) in Eq. (A- 10) with appropriate integration limits gives 


W(0) = j 2 \j/i(s)e ds + J ^\f/ 2 (s)e i 0s ds + J Q i 0s ds + J W 4 ( s ) e ^ 0S ds (A-18) 


Substituting the expressions in Eqs. (A-13)-(A-16) into Eq. (A-18) results in a straightforward 
integration problem. Evaluating the integrals and simplifying, the result is 


W(0) = 


r 6 + fl 2 > 


(3 — 4 cos 0 + cos 20) 


(A- 19) 


The corresponding Taylor series expansion about 0 = 0 is 


W(0) = 


6+e 2 

30 4 


11 


(3- 4 cos 0 + cos 20) ~ 1 0 4 + 

' 7 


720 


23 
15120 


0 6 


(A-20) 


which is identical to Eq. (17). 

For the points near the end of the data record, Eq. (A-l 1) applies with the coordinate s running 
in the negative direction toward the interior of the data, 


Yn-i(0)= f <p N _i(-s)e i 0 ( N s )(-ds) = e i 0N f (Pj(s)ej 0s ds = e j0T / At y*(0) (A-21) 

J 00 J — 00 

where 

<PM=<p N .i(s) (A-22) 
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since the interpolations at the start and end of the record are mirror images of one another. 

Eq. (A-21) shows that computing endpoint corrections for the beginning of the sampled data will 
also provide corrections for points at the end. 

In Eq. (A- 12), W(0) is applied to the discrete points near the beginning and end of the record 
as if they were interior points. The endpoint correction terms make adjustments relative to this to 
give the appropriate calculation for the endpoints. For example, the fourth point of the data record 
(jc 3 ) is involved in the same four interpolations that an interior point is involved in, but also is 

involved in the interpolation in the first interval, (/ 0 , /|), where three points must be used on the 
right side of the interval (jq, x 2 , and x 3 ), along with the initial point x 0 . Figure A-2 illustrates the 

situation. The correction factor for the fourth point therefore represents the contribution of the 
fourth point to the interpolation for the first interval. The Lagrange polynomial multiplying the 
fourth point for interpolation in the first interval is 


<h(s) 


s(.s-l)(s — 2) 

_ 


(A-23) 


Using the appropriate limits, the endpoint correction can be evaluated from 


r ^ 6) = t jdSds = £ 


j($ — 1)($ — 2) 
6 


e ds 


(A-24) 


Evaluating the integral and simplifying, 


r 3 (0) = 


2(3- 0 2 ) - (6 + e 2 )cos 6 
— 


60 — (6 + 0 2 }sin0 
' 60 A 


(A-25) 
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The corresponding Taylor series expansion about 9 = 0 is 


2(3-6> 2 )-(6 + 0 2 )cose 
m ) 604 


60 - (6 + 9 2 }sin0 


60 4 


24 


— e 2 +—^—e A — e 6 -je{— 

180 24192 259200 \360 


(A-26) 

— o 2 + — — — e 4 — — oA 

840 362880 29937600 ) 


which is identical to Eq. (21). 

For y 0 ( 6), Yj ( 0), and y 2 (< &)> the contribution for the first interval interpolation is included in a 
fashion similar to that shown for y 3 (0) above. In addition, for y o (0), Y\(Q)> and y 2 (0), one or 
more terms must be subtracted to account for the fact that x 0 , x x , and x 2 are not involved in all four 
interpolations shown in Figure A-l for an interior point. For example, y 2 (#) includes the 
contribution to the first interval interpolation for x 2 , along with a subtraction for one of the four 
interpolations shown in Figure A-l which does not apply for x 2 . The contribution for this 
interpolation has been included in the W(0) term of Eq. (A-12) for x 2 and must therefore be 
subtracted in the endpoint correction term. The expression for 72 ( 0 ) is then 



s(.y-l)(s-3) 


-2 


e i 0s ds - J Vi(s)e i 0s ds 


(A-27) 


where y/’ 1 (s) is given by Eq. (A-13). Comparing Eqs. (A-27) and (A-ll), 


fc(l)= ^-^-3) 0£J£| 


<Pl W = Vi(») - (^±jX£±jX£±l) -2 <jS-1 

6 


(A-28) 


(p 2 (s) = 0 otherwise 


so that (p 2 (s) is seen to be a function that has non-zero values only for s in the intervals (0,1) or 
(-2,-1). The other endpoint corrections / o (0) and 7,(0) are computed similarly. 
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Figure 7 Finite Fourier Transform Calculation Comparison 
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Figure 10 Finite Fourier Transform Calculation Comparison 
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Fieure 12 Differences for Finite Fourier Transform Calculations 





I I I I 

'i-3 'i-1 'i 


I I I 

t t t * 

i+1 i+2 i+3 


Figure A-l Interior Interpolation 


0 1 


1 o 1 1 


Figure A -2 First Interval Interpolation 


30 


31 



REPORT DOCUMENTATION PAGE 

Form Approved 
OMBNo. 0704-0188 

Public reporting burden for this collection of information is estimated to average 1 hour per response, including the time for reviewing instructions, searching existing data sources, 
gathering and maintaining the data needed, and completing and reviewing the collection of information. Send comments regarding this burden estimate or any other aspect of this 
collection of information, including suggestions for reducing this burden, to Washington Headquarters Services, Directorate for Information Operations and Reports, 1215 Jefferson Davis 
Highway, Suite 1204, Arlington, VA 22202-4302, and to the Office of Management and Budget, Paperwork Reduction Project (0704-0188), Washington, DC 20503. 

1. AGENCY USE ONLY ( Leave blank) 2. REPORT DATE 3. REPORT TYPE AND DATES COVERED 

June 1997 Technical Memorandum 

4. TITLE AND SUBTITLE 

High Accuracy Evaluation of the Finite Fourier Transform Using Sampled 
Data 

5. FUNDING NUMBERS 

522-22-21-01 

6. AUTHOR(S) 

Eugene A. Morelli 

7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES) 

NASA Langley Research Center 
Hampton, VA 23681-0001 

8. PERFORMING ORGANIZATION 
REPORT NUMBER 

9. SPONSORING / MONITORING AGENCY NAME(S) AND ADDRESS(ES) 

National Aeronautics and Space Administration 
Washington, DC 20546-0001 

10. SPONSORING /MONITORING 
AGENCY REPORT NUMBER 

NASA TM- 110340 

11. SUPPLEMENTARY NOTES 

12a. DISTRIBUTION / AVAILABILITY STATEMENT 

Unclassified - Unlimited 
Subject Category - 05 

12b. DISTRIBUTION CODE 

13. ABSTRACT (Maximum 200 words) 

Many system identification and signal processing procedures can be done advantageously in the 
frequency domain. A required preliminary step for this approach is the transformation of sampled time 
domain data into the frequency domain. The analytical tool used for this transformation is the finite 
Fourier transform. Inaccuracy in the transformation can degrade system identification and signal 
processing results. This work presents a method for evaluating the finite Fourier transform using cubic 
interpolation of sampled time domain data for high accuracy, and the chirp z-transform for arbitrary 
frequency resolution. The accuracy of the technique is demonstrated in example cases where the 
transformation can be evaluated analytically. Arbitrary frequency resolution is shown to be important 
for capturing details of the data in the frequency domain. The technique is demonstrated using flight 
test data from a longitudinal maneuver of the F- 18 High Alpha Research Vehicle. 

14. SUBJECT TERMS 

system identification, Fourier transform, frequency domain, accuracy, 
frequency resolution, flight test data 

15. NUMBER OF PAGES 

34 

16. PRICE CODE 

A03 

17. SECURITY CLASSIFICATION 
OF REPORT 

Unclassified 

18. SECURITY CLASSIFICATION 
OF THIS PAGE 

Unclassified 

19. SECURITY CLASSIFICATION 
OF ABSTRACT 

Unclassified 

20. LIMITATION OF ABSTRACT 


NSN 7540-01-280-5500 Standard Form 298 (Rev. 2-89) 

Prescribed by ANSI Std. 239-18 
298-102 


32 


























